#include <sys/time.h>
// Function to get time
double get_cur_time() {
	struct timeval   tv;
	struct timezone  tz;
	double cur_time;

	gettimeofday(&tv, &tz);
	cur_time = tv.tv_sec + tv.tv_usec / 1000000.0;

	return cur_time;
}

__shared__ float3 sharedPos[256];

__constant__ precf D_SIGMA_OO=2.4916;
__constant__ precf D_SIGMA_CC=3.3445;
__constant__ precf D_SIGMA_OC=1.8778;

__constant__ precf D_EPS_OO=7.3638;
__constant__ precf D_EPS_CC=0.85224;
__constant__ precf D_EPS_OC=38.808;

__device__ precf computeDistances(float3 ri, float3 rj, float3 box) { 

	precf r, cp, cx, cy, cz; 

	cx = ri.x - rj.x;
	cy = ri.y - rj.y;
	cz = ri.z - rj.z; 

	//apply periodic boundary conditions
	cx = cx - box.x*rintf(cx/box.x);
	cy = cy - box.y*rintf(cy/box.y);
	cz = cz - box.z*rintf(cz/box.z);

	cp = cx*cx + cy*cy + cz*cz;

	r = sqrtf(cp);

	return r; 

}


__host__ prec radius(int i, int j, int k, int l){

	prec r, cp, cx, cy, cz;

    //get r
	cx = (hAtom[i*3+k].x-hAtom[j*3+l].x);
	cy = (hAtom[i*3+k].y-hAtom[j*3+l].y);
	cz = (hAtom[i*3+k].z-hAtom[j*3+l].z);

	//apply periodic boundary conditions
	cx = cx -box.x*round(cx/box.x);
	cy = cy -box.y*round(cy/box.y);
	cz = cz -box.z*round(cz/box.z);

	cp=cx*cx+cy*cy+cz*cz;
	r = sqrt(cp);

	return r;

}

__host__ prec cpuEnergy(double r, int i, int j, int k, int l){

	prec sig, eps;
	prec dP = 0.0;
	prec a = 2.4829;//cutoff term
	prec B = 2.7956455;

	if(k != 2 && l == 2 || k == 2 && l != 2) {
		sig = SIGMA_OC; 
		eps = EPS_OC;
    }
	if(k != 2 && l != 2) {
		sig = SIGMA_OO;
		eps = EPS_OO; 
	}
	if(k == 2 && l == 2) {
		sig = SIGMA_CC; 
		eps = EPS_CC; 
	}
	prec rcut = a*sig;

	if(r < rcut) {
		prec f1 = sig/r;
		prec f2 = f1*f1;
		prec f6 = f2*f2*f2;
		prec xp = exp(sig/(r-a*sig));
		dP = eps*((B*f6*f6)-f6)*xp;
	}
	else dP = 0;

	return dP;
}

__device__ prec gpuEnergy(prec r, prec sig, prec eps) { 

	precf dP = 0.0;
	precf a = 2.4829f;//cutoff term
	precf B = 2.7956455f;

	precf RCUT = a*sig;

	if(r < RCUT) {

		precf f1 = sig/r;
		precf f2 = f1*f1;
		precf f6 = f2*f2*f2;
		precf xp = expf(sig/(r-a*sig));
		dP = eps*((B*f6*f6)-f6)*xp;
	}

	else dP = 0;

	return dP;
}

__device__ prec phi0(float3* ri, int j, float3 box, float3* hAtom) {

	precf sig, eps; 
	prec dP=0.0; 
	precf r;	

	for(int k=0; k<=2; k++) {
		for(int l=0; l<=2; l++) {
			if (k != 2 && l == 2 || k == 2 && l != 2) { sig = D_SIGMA_OC; eps = D_EPS_OC; }
			if (k != 2 && l != 2) 		    { sig = D_SIGMA_OO; eps = D_EPS_OO; }
			if (k == 2 && l == 2) 		    { sig = D_SIGMA_CC; eps = D_EPS_CC; }

			float3 rj = sharedPos[j*3+l];	
 
			r = computeDistances(ri[k], rj, box);
			dP += (prec) gpuEnergy(r, sig, eps);

		}
	}	

	return dP;
}

__global__ void CO2Energy(float3 *gAtom, prec* d2body, int numMolecules, float3 box) { 
	
	prec twobodyEnergy = 0.0; 
	
	int gtid = blockIdx.x*blockDim.x+threadIdx.x; 

	float3 ri[3]; 
	ri[0] = gAtom[gtid*3+0];
	ri[1] = gAtom[gtid*3+1];
	ri[2] = gAtom[gtid*3+2];

	for (int i=0, tile=0; i < numMolecules; i+=8, tile++) { 

		int idx = tile*blockDim.x + threadIdx.x;

		sharedPos[threadIdx.x*3+0] = gAtom[idx*3+0];   // Carbon
		sharedPos[threadIdx.x*3+1] = gAtom[idx*3+1]; // Oxygen
		sharedPos[threadIdx.x*3+2] = gAtom[idx*3+2]; // Oxygen

		__syncthreads(); 

		for (int j=0; j < blockDim.x; j++) {
			if (gtid < (tile*blockDim.x+j)) {
				twobodyEnergy += phi0(ri, j, box, gAtom); 
			}
		}	 		
		__syncthreads(); 
	}
	
	d2body[gtid] = twobodyEnergy; 	
}
	
		
__host__ prec phi1(int moveSelect) {
	
	prec P1 = 0.0;
	prec dP,r;

	int i = moveSelect;

	for (int j = 0; j<=(nMolecules -1); j++) {
		if (j != i) {
			for(int k = 0; k<=2; k++) {
				for(int l = 0; l<=2; l++) {
					r = radius(i, j, k, l);
					dP = cpuEnergy(r, i, j, k, l);
					P1 = P1 + dP;
				}
			}
		}
	}
	
	return P1;
}

prec calculatePotential(int ref, int moveSelect, float3 *hAtom) {

	prec gpuEnergy_1 = 0.0;

	for (int i=0; i<nMolecules; i++)  {
		h2body[i] = 0.0;
	}

	// Calculate change in energy on the CPU
	if(ref == 1) {
		pot = phi1(moveSelect);
	} else { 	// Compute total energy on the GPU 

		// HOST --> DEVICE MEMORY TRANSFER 
		double startTime = get_cur_time(); 
		cudaMemcpy(dAtom, hAtom, nMolecules*3*sizeof(float3), cudaMemcpyHostToDevice);
	
		int numThreads = 8; 
		dim3 threadDims(numThreads, 1, 1);
		int blockDim = (nMolecules/numThreads);
		if (nMolecules%numThreads != 0) blockDim++;
		dim3 gridDims(blockDim, 1, 1);

		CO2Energy<<<gridDims, threadDims>>> (dAtom, d2body, nMolecules, box);
		cudaThreadSynchronize();

		// DEVICE --> HOST MEMORY TRANSFER

		cudaMemcpy(h2body, d2body, nMolecules*sizeof(prec), cudaMemcpyDeviceToHost);

		for (int i=0; i<nMolecules; i++) {
			gpuEnergy_1 += h2body[i];
        }
		
		double endTime = get_cur_time(); 
		//printf(" Total time is %g\n", endTime-startTime); 

		pot = gpuEnergy_1; 
	
	}

	return pot*ENERGY_CF;

}


